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Abstract. We present the newly developed stochastic model of the galactic cosmic ray 
(GCR) particles transport in the heliosphere. Mathematically Parker transport equation 
(PTE) describing non-stationary transport of charged particles in the turbulent medium is the 
Eokker-Planck type. It is the second order parabolic time-dependent 4-dimensional (3 spatial 
coordinates and particles energy/rigidity) partial differential equation. It is worth to mention 
that, if we assume the stationary case (df /dt = 0) it remains as the 3-D parabolic type problem 
with respect to the particles rigidity R. If we fix the energy {df/dR = 0) it still remains 
as the 3-D parabolic type problem with respect to time. The proposed method of numerical 
solution is based on the solution of the system of stochastic differential equations (SDEs) being 
equivalent to the Parker’s transport equation. We present the method of deriving from PTE 
the equivalent SDEs in the heliocentric spherical coordinate system for the backward approach. 
The obtained stochastic model of the Eorbush decrease of the GCR intensity is in an agreement 
with the experimental data. The advantages and disadvantages of the forward and the backward 
solution of the PTE are discussed. 


1. Introduction 

Many problems in physics, finance or biology can be represented as models of the diffusive 
transport processes described by the Fokker-Planck type equations (FPE). The difficulty of 
the numerical solution of this type equation increases with the problem dimension. Reason is 
the instability of the numerical schemes like finite-differences and finite-volume in the higher 
dimensions. To ensure the scheme stability and convergence the density of numerical grid must 
be improved, increasing the computational complexity. To overcome this problem the stochastic 
methods can be applied (e.g.[I],12]). In this approach the individual particle motion is described 
as a Markov stochastic process, and the system evolves probabilistically. Accordingly, the FPE 
can be solved by the corresponding stochastic differential equations (SDEs) (e.g. [3]). 

We apply stochastic methodology to model the galactic cosmic rays (GCR) transport in the 
heliosphere. During the propagation through the heliosphere GCR particles are modulated by 
the solar wind and heliospheric magnetic field (HMF). Modulation of the GCR is a result of 
action of four main processes: convection by the solar wind, diffusion on irregularities of HMF, 
particles drifts in the non-uniform magnetic field and adiabatic cooling (e.g. H). Transport of 
the GCR particles in heliosphere is described by the Parker transport equation (PTE) [5] being 


the second order parabolic type partial differential equation: 

^ = V • • V/) - (iT, + [/) . V/ + |(V • (1) 

where / = f{f^R^t) is an omnidirectional distribution function of three spatial coordinates 
r — r{r^9^(p)^ particles rigidity R and time t; U is solar wind velocity, the drift velocity, and 
Kf- is the symmetric part of the diffusion tensor of the GCR particles. 

Based on the stochastic approach we model the short time variation of the GCR intensity, 
called the Forbush decrease (Fd) [H]. The Fd occurs as an occasional decrease in GCR intensity 
recorded on the Earth surface by the neuron monitors.The Fds follow the occurrence of the solar 
flares and intensive solar coronal mass ejecta (CME) [TU]. There can be distinguished two types 
of the Fds: 1) sporadic- being the result of the shock waves and the magnetic clouds appearing 
in the interplanetary space, as the result of the solar flares on the Sun; 2) recurrent type Fd 
connected with the corotating interaction regions (CIR) appearing in the interplanetary space 
in connection with the solar rotation. 

2. Stochastic approach 

To model the GCR transport by the stochastic methods the corresponding SDEs must be derived. 
Firstly, the PTE (Eq{T]) must be transformed to the standard form of the FPE which, depending 
on integration’s direction, can be generally expressed in two forms [3]: 
time-for ward: 
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Corresponding to Eqs. [2]and|3]SDE can be written as (e.g. [5]): 

dr = Ai ■ dt + Bij ■ dW, 
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where r is the trajectory of individual pseudoparticle in the phase space and dWi is the Wiener 
process, usually written as dWi = • dwi^ dwi is the randomly fluctuating term having 

Gaussian distribution. 

In both, forward and backward cases, first we start from some initial position in space and time 
and integrate along the pseudoparticles trajectories until they reach the boundary. The choice 
between the forward and backward integration depends on the problem that has to be solved. 
In the case of the GCR propagation in the forward approach particles should be initialized 
at various points at the boundary, where the GCR particles enter the heliosphere. Then, its 
trajectory should be traced until they arrive the point of interest, e.g. Earth orbit. To obtain the 
reasonable statistic a huge number of particles should be initialized because most of them do not 
reach the Earth orbit. The backward approach is more efficient for the GCR propagation in the 
heliosphere, because it reduces the number of ’useless’ particles. In the backward integration 
particles are initialized at the Earth orbit and traced backward in time until they reach the 
heliosphere’s boundary (in this paper assumed at 100 AU, Fig. [T]). The particle distribution 
function can be obtained by averaging over the entrance points, /(r, i?) = ■^Yln=ifLis{R)^ 
where fLis{R) is the cosmic ray spectrum at the outer boundary taken as in [T3j and R is the 




rigidity of the particle at the entrance point. 

The PTE (EqJT|) in the 3-D spherical coordinate system (r, 9, ip) can be written as time-backward 
EPE diffusion equation: 
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with following coefficients: 
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We apply the full 3D anisotropic diffusion tensor of GCR particles Kij = ^ + K^- ’ consisting 




of the symmetric and antisymmetric parts presented in [ 7 ]. The drift velocity of 


GCR particles is implemented as: — 
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[8j. The corresponding to Eq. [5]set of SDEs with 


matrix (i, j = r, 0, (^) has a form (the same form can be found in [E]): 
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The Eqs. [6] are integrated backward in time by the Euler—Maruyama scheme. As an initial 
condition, an empty heliosphere is assumed, as discussed in US]. Solving Eqs. [6] in spherical 
coordinates the following boundary conditions are assumed: (pi < t) ^ ipi = ipi + 27r, 
Pi > 271 ^ Pi = Pi — 277, 9i < 0 ^ 9i = 9i -\- 271 and 9i > 71 ^ 9i = tt — \9i\. For the inner radial 
boundary we assume the reffecting boundary, ^ = 0 at r = 0.001 AU, with time from 0 to 100 
days. The outher boundary is specified at radial distance lOOAU as /(100,i?, t) = fLis{R)- 


3. Model of the short time variation of the GCR intensity 

In this paper we present the model of the recurrent Ed taking place due to established corotating 
heliolongitudinal disturbances in the interplanetary space. CIR passing the Earth gradually 
diminishes the diffusion at the Earth orbit, causing larger scattering of the GCR particles, and in 
effect fewer GCR particles reach the Earth. We simulate this process by the gradual decrease and 
then the increase of the diffusion coefficient at the Earth orbit with respect the heliolongitude. 
The diffusion coefficient Kn of of cosmic ray particles has a form: Kn — Kg • K{r) • K{R^v)^ 
where — lO^^cm^/s, K{r) = 1 + 0.5 • {r/lAU) and K{R^iy) = R^~^ . The exponent u 
pronounce the increase of the HMF turbulence in the vicinity of space where the Ed is created 
(e.g. [lU [E]), and is taken as: u = 0.8 + 0.2sin{p — 90°) for 90° < p < 270°. We assume 
in the model the existence of the two dimensional spiral Parker’s heliospheric magnetic field 
m implemented through the angle ^ = arctan{—By:^/Br) = arctan{yt • r • sin9jU) in the 3D 
anisotropic diffusion tensor Kij of GCR particles [7]. 

The expected changes of the GCR intensity for the rigidity of 10 and 20 GV during the simulated 
Ed in comparison with the profiles of the daily GCR intensities recorded by the two neutron 


































Figure 1. The sample pseudoparticles 
trajectories within the heliosphere. 



Figure 2. Changes of the expected 
amplitudes of the Fd of the GCR intensity 
at the Earth orbit, for the rigidity of 10 and 
20 GV based on the solutions of the backward 
SDEs in comparison with the GCR intensity 
registered by Moscow and Beijing neutron 
monitors during the Fd in March 2002. 


monitors with different cut off rigidities in 18 March - 4 April 2002 presents Fig. [21 One can 
see that as is expected the amplitude of the Fd decreases for higher rigidities. One can see that 
the proposed model is in a good coincidence with the experimental data. Moreover, the model 
of the Fd obtained based on the solution of the SDE allows to reffect the stochastic character 
of the GCR particles distribution in the heliosphere and present the pseudoparticle trajectory 
thorough the 3D heliosphere (Fig. [T]), which is not possible based on the solution of the Parker 
transport equation by the e.g. finite difference method (e.g. m)- 
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